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ABSTRACT 

In the galaxy clusters SSDS J1004+411 and AGO 1689, a large number of 
multiply-imaged background objects have recently been observed. We use these 
data to map the projected mass distribution in the inner regions of these clusters. 
The source redshifts in these clusters span a large range, which eliminates the 
degeneracies that plague nearly all lensing work. As a result the mass maps are 
very well-constrained, despite very weak prior assumptions. AGO 1689 lenses 
so many objects that we can afford to map it twice using disjoint sets of im- 
ages, thus verifying our internal error estimates. We then deproject the mass 
maps (pretending for this purpose that they are spherical) and obtain inner pro- 
files consistent with p oc r~^ and indistinguishable from recent cold dark matter 
simulations. Assuming that baryons make only a small difference to the profile 
outside ~ 10 kpc, these results support the prediction of a universal inner profile 
for dark matter structures, independent of any parametrization of that profile. 

Subject headings: gravitational lensing; clusters: individual (SDSS J1004+411, 
AGO 1689); dark matter 



1. Why inner profiles of clusters? 

A key prediction of heirarchical clustering dominated by cold dark matter is that viri- 
alized structures have a universal density profile. Although there is still no first-principles 
derivation of a universal profile, A^-body simulations and phenomenological models both 
consistently indicate that such a profile exists. 
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The best-known candidate for a universal profile remains the original proposal 



p oc (r/r,) ^(1 + r/r,) 



(1) 



from Navarro et al. (1997), but it is gradually being superseded. In recent simulations — see 
for example, Fig. 3 of Navarro et al. (2004) — NFW remains a reasonable fit to individual 
simulated halos, but an ensemble of halos shows the steepness changing gradually with no 
distinct r^. Graham et al. (2005) discuss alternative parametrizations. The smallest resolved 
scales in current dark-matter simulations have p ~ j--i-2±o.2 (^Qiemand et al. 2004). The 
total mass profile would be somewhat steeper at the smallest radii, because of the infall of 
baryons and the response of the dark-matter to the baryons. The latter effect is generally 
modelled as adiabatic contraction (Blumenthal et al. 1986; Gnedin et al. 2004). 

These results motivate research on the inner profiles of real clusters. 



Testing predicted universal profiles against observations seems straightforward at first. 
Well-known strategies are to 

1. put observed velocity dispersions through the Jeans equation to infer a density profile. 



2. fit a mass model to lensing data, or 

3. invert an observed x-ray profile through the hydrostatic equilibrium equation. 

The first two approaches, and to a lesser extent the third, have all been followed in the 
literature, but the conclusions remain unclear. For example: Lokas et al. (2006) using 
method (1) find NFW profiles satisfactory, whereas Sand et al. (2004) using a combination 
of (1) and (2) find inner profiles significantly shallower than NFW. 

Evidently there are major uncertainties involved in extracting a mass profile. It is useful 
to get some insight into what the uncertainties are, and hence how robust different methods 
are likely to be. By considering the basic equations involved, we can identify problematic 
operations that either amplify noise or introduce systematic errors. Let us call them gremlins. 

The first gremlin appears in deprojecting an observed quantity, via the well-known Abel 
integral 



2. Why J1004+411 and AGO 1689? 



or 




(2) 
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In this case, differentiating an observed amplifies the noise. Smoothing will 

suppress the noise of course, but will replace it with a systematic error. However, the 
integral tends to suppress the noise somewhat, so this gremlin is comparatively benign. 

Next, consider Jeans equation for spherical systems 

GM{r)^'j^^{p,a'^)-2r{a'^-a^) (3) 

relating the luminous density pi and velocity dispersions to the mass profile. Since cr^ and ctj 
cannot be directly observed, they need to be related to the projected dispersion somehow. 
This introduces a second gremlin: assumptions about anisotropy produce a systematic error; 
deriving the anisotropy from higher velocity moments (Lokas et al. 2006) introduces noise. 
The derivative in Eq. (3) is a third gremlin. 

For x-ray emitting gas in hydrostatic equilibrium, the temperature profile T(r) and gas 
density profile Pg(r) are related to the total mass profile by 

GM(r) = -— + P 4) 

/inip \T dr Pg dr J 

where Pp is the proton mass and p the particle weight. Here the first gremlin appears in 
deprojecting for T{r) and Pg(r), while the third gremlin appears in another guise in the 
derivatives of T and Pg. 

Lensing is free of the second and third gremlins. Moreover it provides one simple and 
robust constraint on mass profiles. Central images are extinguished by an isothermal or 
steeper cusp (Evans & Wilkinson 1998), and the fact that clusters commonly host central 
images implies that clusters are shallower than isothermal (Subramanian et al. 2000). 

Lensing, however, introduces a new gremlin, namely lensing degeneracies. Consider the 
equation for the arrival time^ 

r(0,; 0s) = \{e,- e,f - 2c,sV-2«;(0O (5) 

which represents the time delay for a virtual photon^ originating from a source at sky position 
Os and getting deflected at the lens so that the observer sees it coming from sky-position 



^Eq. (5) is Eq. (2.6) from Blandford & Narayan (1986), except that we have introduced three new symbols: 
V^^ is an operator that solves Poisson's equation in 2D, Cls — dj^s/d-os, while k = ATrdoj^Y] is the surface 
density in units of the critical density for sources at infinity. Note that Blandford & Narayan set G = c = 1. 

^Real photons take paths for which r is a maximum, minimum, or saddle point (Fermat's principle). 
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Oi. The first term on the right of Eq. (5) is the geometrical path difference between a 
deflected and undeflected photon trajectory, and the last term is the Shapiro time delay. 
The whole equation is dimensionless: r is a scaled time delay and k is a scaled surface 
density, the coordinates are in radians, while 

~ 1 is a ratio of distances depending on 
lens and source redshifts. The gradient of r with respect to 0i gives the lens equation, 
and the matrix of second derivatives is the inverse magnification. Thus, image positions 
in strong lensing constrain the first derivative of r while shear or magnification in weak or 
strong lensing constrain the second derivative. These derivatives are not gremlins, because 
no differentiation of data is involved. But let us rewrite the equation as 

T = 2v-\i - c^s^i) - ■ ds (6) 

where we have expanded the geometrical part and then discarded a 0g term since it has no 
optical effect. If Cls is a constant (that is, there is only one source redshift) then we may 
rescale Eq. (6) by some positive factor A: 

T ^ Ar, ds^XOs, (1 - Clsk) ^ A(l - Clsk). (7) 

Eq. (7) amounts to redefining the contour spacing in Blandford & Narayan's Fig. 2, without 
changing the figure — the observables don't change, but the k. profile becomes steeper (if 
A > 1 or shallower (if A < 1). This is our fourth, and perhaps most vicious, gremlin. The 
case where A is constant is well known as a steepness degeneracy (Falco et al. 1985), but 
in fact A need not even be constant (Saha 2000; Zhao & Qin 2003; Saha & Williams 2006). 
Having a range of source redshifts breaks the degeneracies (Abdelsalam et al. 1998b); having 
number counts of weakly lensed objects breaks the main steepness degeneracy (Dye et al. 
2002). 

The above arguments indicate that lens reconstructions with special attention to break- 
ing lensing degeneracies are hkely to provide the most robust mass profiles. The clusters 
which have been so analyzed are AGO 2218 (Kneib et al. 1996; Abdelsalam et al. 1998b) 
AGO 370 (Abdelsalam et al. 1998a), AGO 1689 (Dye et al. 2002; Broadhurst et al. 2005; 
Diego et al. 2005; Zekser et al. 2006; Halkola et al. 2006), and J1004+411 (Sharon et al. 2005). 
The first two have double-peaked mass distributions, evidently from ongoing mergers. That 
leaves the last two lensing clusters as worthy of special attention. 

3. Mass profiles 

We will now use the excellent recent strong-lensing observations by Sharon et al. (2005) 
on J1004-I-411 and by Broadhurst et al. (2005) on AGO 1689 to map the mass distribution 
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in the inner regions of these clusters. As we noted above, these data have aheady made 
into models by the observing authors and others. Our work, however, goes further than 
previous models because (i) rather than fitting one or a few models to the data, we gener- 
ate large ensembles of models exploring the possible mass distributions that can reproduce 
the data; (ii) we compare results from independent data sets to test for consistency; and 
(iii) we deproject the surface density profiles non-parametrically for a better comparison 
with simulations. 

Our technique is implemented in the PixeLens code, originally developed for galaxy 
lenses (Saha & Williams 2004), but now enhanced to use multiple sourcc-rcdshifts. Given 
the multiple-image data, the program generates an ensemble of 400 models that each (a) re- 
produce the lensing data exactly, and (b) are consistent with a given prior on the mass map. 
The model ensemble automatically provides Bayesian uncertainties on the mass map or any 
derived quantity. We assume a concordance cosmology {K = 0, Jl^ = 0.3, ty = — 1) with 
Hq^ = 14 Gyr. Cosmological parameters enter into the the mass and distance scales, but 
not the profiles themselves. 

The prior we assume is the same as in Saha & Williams (2004) — basically it requires 
that that the surface density be non-negative and centrally concentrated, with substructure 
allowed — but with one difference: we remove the requirement that E be steeper than 
(which was motivated by galaxy dynamics) and require only that the cylindrically-averaged 
S decrease with R. We measure R from the center of the brightest cluster galaxy, so the 
mass peak is modeled to coincide with the light peak. In some lensing clusters mass and 
light peaks have been inferred to be offset tens of kpc (Abdelsalam et al. 1998a,b), but in 
J1004 -1-411 and AGO 1689 the central lensed images indicate that any such offset is small. 

The requirement of exact fits to the image data leaves PixeLens currently unable to 
handle all the multiple-image data on AGO 1689. (The problem is accumulation of roundoff 
error.) Nevertheless, we can make a virtue of this by splitting the data into two independent 
sets. Accordingly, we present three mass reconstructions: 

1. J1004+411 (zj^ — 0.68) reconstructed using 13 images coming from 4 sources (zg = 1.73 
to 3.33), 

2. AGO 1689 (zl = .18) reconstructed using 30 images of 9 sources {zs — 1.74 to 4.52), 
and 

3. AGO 1689 again, reconstructed using 28 images of 7 other sources {zs = 1.57 to 5.16). 

Figure 1 packs all the results of this paper. To the left we show three ensemble- average 
mass maps. Both clusters appear fairly round, but have some substructure, such as a fea- 
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ture to the NE in AGO 1689. In this paper we concentrate on radial dependence, leaving 
substructure for later work. 

The middle column of the figure shows the cylindrically- averaged surface density 
with 90%-confidence uncertainties derived from the model-ensemble. In each case we see a 
shallow and tightly constrained profile in the inner region, and moreover the two profiles for 

AGO 1689 from independent image-sets agree. (The profile beyond the outermost image is 
to be disregarded; there PixeLens has no data to work with and simply makes the profile fall 
smoothly to zero at the edge of the map.) The enclosed mass agrees with previous strong- 
lensing estimates but is about a factor of two higher than current models for the x-ray gas 
give (Andersson & Madejski 2004). The improvement upon our earlier reconstructions of 
J1004 -1-411 WiUiams & Saha (2004) using a single 4-image background object is remarkable: 
then we concluded that the lens was consistent with NFW, but even a profile as steep as 
R"^-^ was not excluded; now simply inspecting the upper middle panel of Fig. 1 shows that 
this cluster is shallower than 1/R. 

To the right of the figure we show p(r), derived from S(i?) using Eq. (2). To evaluate the 
numerical derivative and then the integral, we fit a spline to E up to the R of the outermost 
image and assumed E oc outside. (The result is not sensitive to the assumed outside 
slope, or to the boundary between spline and provided the latter is not beyond the image 
region.) Here again we show ensemble-derived uncertainties, which in this case are probably 
too small, because substructure will introduce systematic errors through the derivative in 
Eq. (2). With this caveat, our constraints on the steepness between 10 kpc and the outermost 
image are as follows: J1004-I-411 is best fit by p oc r■^^^^ but the uncertainty bands shown 
allow r~°-^ to r~^-^, while AGO 1689 is best fit by p oc r~^-'^^ but the uncertainty bands allow 
^-1.0 j,-i.25 xhese values are typical of the inner slopes in Table 3 of Diemand et al. (2004). 
The right column of Fig. 1 is also comparable to their Fig. 4. For AGO 1689 we have also 
attempted to model the change in the profile due to adiabatic contraction, using estimated 
gas density from Fig. 9 of Andersson & Madejski (2004), but the difference is smaller than 
our estimated uncertainties. It appears that baryons make a difference of at most 10%, which 
is of the same order as the substructure. We conclude that, within current uncertainties, the 
inner slopes in these lensing clusters and in simulations are indistinguishable. 



4. Discussion 

In this paper we have used recent strong lensing observations having multiple source- 
redshifts (Sharon et al. 2005; Broadhurst et al. 2005) to reconstruct the inner density profiles 
of two clusters in a non-parametric way. These inner density profiles show remarkable agree- 
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ment with simulations of CDM halos. Although we have disregarded the secondary questions 
of baryons (not included in simulations at this scale) and substructure (present in both the 
real Universe and simulated ones), our results arguably confirm an important prediction of 
heirarchical structure formation. 

But why are the profiles the way they are? There is still nothing approaching a first- 
principles derivation of a universal profile, though there arc several interesting ideas. The 
only exactly solved case is self-similar spherical collapse, but that implies p oc r~^'^^ (Fillmore 
& Goldreich 1984; Bertschinger 1985) which is too steep. Subramanian et al. (2000) propose 
a cunning variant involving merging self-similar proto-halos into a Matreshka-likc hcirarchy, 
resulting in a composite halo with a shallow inner slope; but their mechanism relies on undi- 
gested merging. Several other scenarios involving spherical collapse have also been studied, 
and recently Lu et al. (2006) proposed a phenomenological ID model that produces an 
inner slope and an outer slope. Finally, one numerical result is especially intriguing: in 
CDM collapse simulations the surrogate phase-space density p/a^ oc r~° (Taylor & Navarro 
2001; Barnes et al. 2005). Combining this relation with the Jeans equation gives p oc r~°'* 
asymptotically for small r (Austin et al. 2005; Hansen et al. 2006; Dehnen & McLaughlin 
2005). The total mass profile at very small r would be steeper than this because of the 
baryons, and anyway our results do not have the resolution to test for this yet. But it 
remains an interesting prediction. 
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Fig. 1. — Mass maps and profiles for J1004+411 and for AGO 1689. Left column: Ensemble- 
average mass maps together with the lensed image positions. Panels are oriented as North- 
up/East-left. Ticks correspond to 20 kpc for J1004+411 and 100 kpc for AGO 1689. Gontour 
steps are 0.5 x lO^M© kpc^^ for J1004+411 and IO^Mq kpc^^ for AGO 1689. Middle 
column: for the mass maps in in the left column, together with 90% uncertainties. 

The vertical lines to the right show R for the outermost image. Right column: Density 
profiles, derived by deprojecting the profiles from the middle column. The sloped line shows 

— 1 2 

p oc r . 



